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Using the Crank-Nicholson method, we study the evolution of a Bose-Einstein condensate in an 
optical lattice and harmonic trap. The condensate is excited by displacing it from the center of the 
harmonic trap. The mean field plays an important role in the Bloch-like oscillations that occur after 
sufficiently large initial displacement. We find that a moderate mean field significantly suppresses 
the dispersion of the condensate in momentum space. When the mean field becomes large, soliton 
and vortex structures appear in the condensate wavefunction. 
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I. INTRODUCTION 

The dynamics of Bose-Einstein condensates (BECs) in 
optical lattices have been studied using both theoretical 
and experimental methods (for a review, see Ref. [1]). In 
many cases, the atomic interaction cannot be ignored due 
to the combination of low temperature and high density 
found in a BEC. This interaction is often accounted for 
as a mean- field effect. The role of the mean field in opti- 
cal lattices has received theoretical attention in a number 
of phenomena, such as pulsed atom lasing in an optical 
lattice 0,0]: number squeezing of a BEC in an optical lat- 
tice Q, Landau-Zener tunneling of a BEC [5, 6], and the 
dynamics of a BEC in an accelerating optical lattice 0]. 
In related experimental work, mean-field-induced four- 
wave mixing of matter waves has been observed [8], and 
the mean field has been shown to change the Landau- 
Zener tunneling rate @, 0]. In this paper, we simulate 
the dynamics of a BEC in a combined potential con- 
sisting of an optical lattice and a harmonic trap. The 
initial displacement of the BEC from the center of the 
harmonic trap triggers motion similar to a Bloch oscilla- 
tion. Bloch oscillations of BECs in optical lattices have 
been studied theoretically M,[Io[ [Tlj and observed exper- 
imentally [HEl QUE 

Il6lll7| in various configurations. 
In this paper, we focus on the systematic investigation of 
the role of the mean field in Bloch oscillations in harmonic 
potentials. We find that the stability of the BEC during 
the oscillations depends on the interplay between the in- 
homogeneity of the applied force and the mean field. In 
recent experiments on Bloch oscillations driven by a con- 
stant external force, the mean field was shown to cause 
momentum dispersion of the BEC [HI, [l6[ • Interestingly, 
in Bloch oscillations in an approximately harmonic mag- 
netic trap, this dispersion can be suppressed by the mean 
field [17]. All of these findings are consistently described 
by the results in this paper. 

The paper is organized as follows. In Sec. [Ill we dis- 
cuss the methods and formalism used in this work. In 
Sec. HIH we explore the effect of the mean field on dis- 
persion characteristics in an inhomogeneous force field. 
In Sec. HV] we discuss solitions and vortices observed in 



the case of a large mean field. In Sec. [Vj we conclude the 
paper. 



II. INTEGRATING THE GROSS-PITAEVSKII 
EQUATION 

At zero temperature and in the limit of weak excita- 
tion, the BEC wavefunction, \£(r, £), can be described by 
the Gross-Pitaevskii (GP) equation [l8|, Qjjj], 



ih 



dt 



2M 



V 2 r + V ext (r)+NV int \*(r,t)\ 2 ]*(r,t) , 

(1) 

where M is the atomic mass, N is the total atom number, 
V ext (r) is the external potential and T^ nt characterizes 
the strength of the mean-field interaction between the 
atoms, defined as Vj nt = 47rft 2 a/M, with a being the s- 
wave scattering length of the atom. 

We assume that the BEC is confined in a cylindri- 
cally symmetric harmonic trap, which can be written as 
Hrap( r ) = \M(u 2 r 2j ruo 2 z 2 ) using cylindrical coordinates 
r = (r, z). In addition, an optical lattice is applied 
to the BEC in the direction of the axis of symmetry, 
VoL(r) = sin 2 (27rz/A), where A is the wavelength of 
the lattice laser. In the case of zero angular momen- 
tum about the symmetry axis, the GP equation can be 
rewritten in cylindrical coordinates as 



dV(r,z,t) _r H 2 1 d d_ _cf_ 
1 dt "^2MV^ r * + &2 ) 

+ iMKV+c, z V) + F sin 2 (^) 
+ NV int \y(r,z,t)\ 2 }V(r,z,t). (2) 

The normalization condition is given by 
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(3) 



In this paper, we numerically solve this 2-dimensional 
(2D) GP equation to study the effect of the mean field, 
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NVi n t\^(r, z,t)\ 2 , on the dynamics of the BEC in the 
optical lattice. 

A well-established numerical method for solving the 
GP equation is the Crank-Nicholson algorithm [20[. To 
arrive at boundary conditions that are conducive to solv- 
ing the problem using the Crank-Nicholson method, we 
introduce a wavefunction y?(r, z,t), defined as 



<p(r,z,t) 



The boundary conditions for <p(r,z, t) are 
(p(r,z,i) — > when r — > or oo, or when \z\ - 
The GP equation in term of (f(r,z, t) becomes 
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that 
» oo. 
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<p(r,z,t) . 



(5) 



To reduce the 2D GP equation into tractable ID steps, 
it is split into separate radial and axial equations: 
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dtpjr, z, t) 
dt 
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<p(r,z,t) 



(6) 



ih 



d(p (r,z,t) 
dt 



h? d 2 1^02 t T • 2 / x 

-wa? + 2 M ^ a + W(— ) 



2 r 



<p(r,z,t) . 



(7) 



In each time increment of the integration, the evolution 
of <p(r,z, t) is broken up into two steps. In the first step, 
we fix the z coordinate and solve Eq. [6] using the stan- 
dard ID Crank-Nicholson method. In the second step, 
the resultant (p(r, z, t) is entered into Eq. [71 as the ini- 
tial wavefunction to calculate the time evolution in the 
^-direction while the r coordinate is fixed, again using 
the ID Crank-Nicholson method. The time and spatial 
increments used in our simulations are 100 ns and 50 nm, 
respectively. A detailed discussion of this variation of the 
Crank-Nicholson method for solving the 2D GP equation 
can be found in Refs. [H [H, Q . 



III. DISPERSION OF THE BEC DURING 
BLOCH OSCILLATIONS 

In close analogy with our recent experiment [l7j , we 
consider a BEC with a fixed number of 87 Rb atoms con- 
fined in a cylindrically symmetric harmonic trap, with 
trap frequencies uj v = 2tt x 20 Hz and co z = 2tt x 50 Hz in 
the radial and axial directions, respectively. In the sim- 
ulations, the atoms are adiabatically loaded into an op- 
tical lattice with a depth of one recoil energy, formed by 
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FIG. 1: (a) Initial BEC wavefunction for 25000 atoms (linear 
gray-scale representation), (b) Relative position of the BEC 
with respect to the harmonic trap along the z-direction at 
t=0. 



two counter-propagating laser beams in the z-direction 
with wavelength A= 852 nm. Then, at time t = 0, the 
harmonic trap is suddenly shifted in the z-direction. To 
obtain the wavefunction immediately before the shift, we 
start with the BEC wavefunction in the 3D harmonic po- 
tential without optical lattice in the absence of a mean 
field, given by the well-known ground state of the har- 
monic potential. The ground state of the BEC including 
optical lattice and mean field is then obtained by grad- 
ually ramping up the optical-lattice and mean-field po- 
tentials to their full strengths over an interval of 100 ms, 
using the formalism described in Eqs. (|6|) and ([7j). The 
initial BEC wavefunction for 25000 atoms is shown in 
Fig. □(a). 

After obtaining the initial wavefunction, the magnetic 
trap is shifted at time t = in the z-direction by a 
distance of Az, leading to the initial condition shown 
in Fig. [U (b). The resulting BEC motion depends on 
Az, and can be broken into several regimes. If Az is 
small enough that the BEC never gains enough momen- 
tum to come close to the edge of the first Brillouin zone, 
the BEC performs small oscillations about the minimum 
of the harmonic trap. For large Az, the BEC peri- 
odically traverses the entire first Brillouin zone. The 
BEC motion is similar to a Bloch oscillation, except 
that the BEC experiences an inhomogeneous rather than 
a constant force. The critical displacement, Az cr , for 
Bragg reflection and Bloch oscillations to occur follows 
from the requirement that the BEC must reach at least 
one recoil energy, £?r, when passing through the min- 
imum of the harmonic trap (with no lattice). Under 
our conditions, E R = h 2 /(2MX 2 ) = h x 3.16 kHz and 
Az cr = ^/2E^/ (Muo 2 ) = 17.1 /am. There is also an in- 
termediate regime in which only a portion of the BEC 
undergoes Bragg reflection. In this paper, we focus on 
the dynamics of the BEC when Az > Az cr and study 
how the mean field affects the BEC during Bloch oscilla- 
tions. 

The initial displacement Az is set to be 30 /am. To 
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FIG. 2: Momentum distributions of the BEC as a function 
of time for different atom numbers, TV. (a)-(d) correspond to 
TV=0, 2000, 25000, and 60000, respectively. The brightness 
represents the atomic density on a linear gray-scale. 
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FIG. 3: ^-momentum probability distribution at t=0 and 
t=15 ms for TV=0, 2000, 25000, and 60000, respectively. The 
center of the p z distribution at t=15 ms is shifted to p z =0. 
The values of a are given in units of p re c- 



visualize the BEC dynamics, we perform a Fourier trans- 
form [24] on the wavefunction \P(r, z, t) and plot the mo- 
mentum distribution of the BEC as a function of time. 
Four different atom numbers are used in the simulations. 
The results are shown in Fig. [2] (a)- (d) for TV=0 (exper- 
imentally, this would be equivalent to a vanishing scat- 
tering length), 2000, 25000, and 60000, respectively. The 
radial (axial) momentum spread of the BEC is given by 
the vertical (horizontal) width of the spots in Fig. [2j The 
momentum plots at different times t are stacked in the 
vertical direction. In the following discussion, we focus on 
the z-degree of freedom, in which most of the interesting 
dynamics occur. As seen in Fig. El the BEC periodically 
scans the first Brillouin zone and is Bragg-reflected at 
times when it reaches one recoil momentum, p rec . Since 
the BEC experiences a near-constant force during the os- 
cillations, its momentum varies at a near-constant rate 
and the overall motion resembles that of a Bloch oscilla- 
tion. 

Comparing the Bloch oscillations for different atom 
numbers, shown in Fig. [2] (a)- (d), we notice that without 
the mean field (TV=0 in Fig. El (a)), the BEC displays a 
significantly larger amount of momentum dispersion than 
in the cases when the mean field is included. While the 
mean field with TV = 25000, shown in Fig.[2](c), seems to 
greatly suppress this momentum dispersion, mean-field- 
induced structure starts to appear in momentum space 
for TV = 60000 in Fig. EI (d). 

To better view the momentum dispersion, we plot the 
atomic density as a function of p z after integrating over 
p r , the momentum in the r-direction. The results are 



shown in Fig. [3] (a)-(d), corresponding to t=0 ms and 
15 ms in Fig. El (a)- (d), respectively. We evaluate the full 
width half maximum, <j, of each curve using a Gaussian 
fit. In Fig. [3] (a), with TV = 0, a increases by a factor of 
more than five from t=0 ms to t=15 ms. This dispersion 
is due to the fact that the force across the condensate is 
not uniform. Atoms farther away from the center of the 
harmonic trap experience a larger force during the oscil- 
lation, thus reaching p rec earlier than atoms closer to the 
center. This effect initiates a breathing motion of the 
BEC wavepacket in coordinate space and causes the mo- 
mentum of the BEC to spread as well. This explanation 
is supported by an additional simulation in which, for 
TV = 0, the z-component of the harmonic trap is replaced 
by a constant force. In this case, the BEC wavepacket 
does not disperse during Bloch oscillations in either co- 
ordinate or momentum spaces. The latter result is con- 
firmed by recent experiments elsewhere [HI, [l6| , in which 
more than 10000 Bloch oscillations are observed when 
the BEC experiences a constant external force and the 
atomic interaction is tuned to through a Feshbach res- 
onance. 

The momentum dispersion is greatly suppressed when 
a moderate mean field is included in the simulation. For 
TV=2000, shown in Fig. [3] (b), a increases by a factor of 
less than two over 15 ms. When TV reaches 25000, there 
is no significant momentum dispersion over 15 ms, as 
can be seen in Fig. [3] (c). To qualitatively explain these 
results, we consider the mean field potential V m f within 
the Thomas- Fermi approximation. Before the shift, the 
BEC has a spatially constant chemical potential, /i, given 
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FIG. 4: (a) Same plot as Fig. [2] for N = 25000, with a 
constant external force applied to the BEC. (b) The harmonic 
potential in the ^-direction is first shifted by Az = —30 /im, 
and then switched to an anti-trap, (c) In addition to (b), 
the harmonic potential in the r-direction is also switched to 
an anti-trap and the interaction between condensate atoms is 
changed to be attractive. 

by 

V = V m t + ^M(cu 2 z 2 +u; 2 r 2 ). (8) 

Therefore, V m { = fi - ±M(lj 2 z 2 + uo 2 r 2 ) inside the BEC. 
Immediately after the shift, as shown in Fig. [T] (b), the 
BEC experiences a total potential Vtot 

Vtot =li+ ^Mco 2 Az 2 - Mco 2 zAz . (9) 
The total force on the BEC, F to t, is given by 

Ftot = - Wtot = Mco 2 Az z . (10) 

According to Eq. [10J with the mean field included, the 
BEC experiences a constant total force when displaced 
from the center of the harmonic trap. In this way, the 
mean field helps stabilize the BEC wavepacket during 
the Bloch oscillations. The validity of the Thomas- Fermi 
approximation requires the mean-field potential to dom- 
inate the external potential. In our case, the approxima- 
tion is valid for N = 25000 but not for N = 2000. 

Based on the above argument, when a constant exter- 
nal force is applied to the BEC, the mean field will cause 
momentum dispersion. To test this, we simulate Bloch 
oscillations under a constant force of the same magni- 
tude as that caused by a Az = 30 fim shift of the trap. 
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FIG. 5: BEC spatial density at (a) 15 ms after a 30 /mi 
displacement, (b) 14 ms after a 20 /im displacement, and 
(c) 10 ms after a 17.2 /im displacement. The atom number 
iV=60000 in all cases. (Linear gray-scale representation.) 



The corresponding simulation results are shown in Fig. [H 
(a), where momentum dispersion is clearly seen. This 
has also been observed experimentally elsewhere [TBI. fl6j. 
To expand on our argument, we simulate the dynamics 
of the BEC under two other conditions with the same 
initial wavefunction as before. In the first case, at t=0 
the harmonic potential in the z-direction is first shifted 
by Az = —30 jam and then switched to an anti-trap. 
Since the mean field now enhances the inhomogeneity of 
the effective force acting on the BEC, we expect faster 
momentum dispersion, which is verified by simulation re- 
sults presented in Fig. 0] (b). The second case has the 
same initial conditions as the first, except that the har- 
monic potential in the r-direction is also changed to an 
anti-trap and the atomic interaction is made attractive. 
In this case, the momentum dispersion is expected to 
be suppressed again because the total force due to the 
harmonic potential and mean field is again constant. In 
Fig. [H (c), no momentum dispersion is seen over several 
Bloch oscillations, as anticipated. Eventually, the BEC 
oscillation becomes unstable due to the intrinsic instabil- 
ity of BECs with negative scattering length. 
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FIG. 6: Same plot as Fig. for N = 60, 000, with (a) Az = 
20 /im and (b) Az = 17.2 /im. Time interval between adjacent 
images is 1 ms. 



IV. CREATION OF SOLITONS AND VORTICES 

For the conditions studied in this paper, an increase 
of the mean field beyond values corresponding to N ~ 
25000 causes discrete peaks to appear in the BEC mo- 
mentum distribution, as seen in Fig.|2](d). This indicates 
that the internal structure of the BEC wavefunction has 
been disturbed. We plot the corresponding BEC den- 
sity, |^(r, z, t = 15 ms)| 2 , in Fig. [5] (a), and can see dark 
stripes within the BEC. These dark stripes are solitons 
created by the BEC standing wave formed during Bragg 
reflection [11]. Alternatively, considering the structure of 
the lowest energy band of the optical lattice, the soliton 
formation can be attributed to the anomalous disper- 
sion [25] or negative effective mass [lH, H?} of the BEC 
near the band edge. The solitons decay into vortices 
as the oscillation goes on. The creation of solitons and 
vortices during Bloch oscillations was first reported in 
Ref. H3. In this section, we show that the initial dis- 
placement Az has a considerable impact on this phe- 
nomenon. 

Figures [6] (a) and (b) show the BEC momentum dis- 
tributions for Az = 20 /im and Az = 17.2 /im, respec- 
tively, both with TV=60000. As shown in Fig. [6] (a), the 
BEC momentum starts to develop structure along the z- 
direction after a single Bragg reflection. In Fig. [6] (b), a 
single Bragg reflection causes even more structure in the 
BEC momentum distribution. Additionally taking into 



account the results shown in Fig. [51(d), we conclude that 
a smaller Az leads to a faster creation of solitons and vor- 
tices. This observation can be explained as follows. Since 
it is the standing wave formed during the Bragg reflec- 
tion that leads to the creation of solitons and vortices, it 
is reasonable to assume that the standing wave duration, 
AT, facilitates the soliton formation. The duration AT 
can also be thought of as the Bragg reflection time, which 
is determined by how fast the BEC passes through the 
anti-crossing region between the two lowest energy bands 
of the optical lattice [17]. Thus, AT is proportional to 
1/F, where F is the force applied to the BEC. Smaller 
displacements Az (under the condition Az > Az cr ) lead 
to smaller F and larger AT, which causes the solitons 
and vortices to appear earlier, as observed in the simula- 
tions. 
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FIG. 7: (Color online) (a) and (c) are zoomed-in images of 
Fig. [5] (b) and (c), respectively, (b) and (d) are phase plots 
corresponding to (a) and (c). (black=0, white=27r) 



Around 14 ms in Fig. El (a) and 10 ms in Fig.[6](b), the 
solitons start to decay into vortices. The corresponding 
BEC densities in coordinate space are shown in Figs. [5] 
(b) and (c), with zoomed-in images presented in Figs. [3 
(a) and (c), respectively. To verify the vortex formation, 
we plot the phase of the wavefunction in Figs. [71 (b) and 
(d). Two of the vortices in each of Figs. 0(b) and (d) are 
highlighted by circles, where the phase changes from 
to 2tt in one round trip. The arrow represents the direc- 
tion of the vortex. Since the vertical axis corresponds to 
the radial direction, and the system is azimuthally sym- 
metric, the two highlighted vortices are part of a single 
annular vortex ring around r = 0. 
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V. CONCLUSION 

In this paper, we have simulated the dynamics of a 
BEC in a combined potential consisting of a magnetic 
trap and an optical lattice. After a sufficiently large ini- 
tial displacement from the center of the harmonic trap, 
the BEC undergoes an oscillatory motion similar to a 
Bloch oscillation. We find that the mean field plays an 
interesting role during these oscillations. While a mod- 
erate amount of mean field suppresses the momentum 
dispersion of the BEC, at large values it causes solitons 
and vortices to appear. We have also investigated the 
dependence of vortex formation on the initial displace- 
ment Az. Some of the results presented in Sec. IHIIhave 



already been verified in a series of Bloch-oscillation ex- 
periments [TBI, [HI, 03 • Quantitative studies of the mean- 
field-induced dispersion characteristics of Bloch oscilla- 
tions in harmonic traps may be performed by time-of- 
flight measurements of the momentum distribution vs 
atom number. It will be challenging to experimentally 
probe the formation of annular vortices in 2D Cartesian 
projections of the BEC obtained using standard time-of- 
flight shadow imaging. However, it may be possible to 
obtain indirect evidence, such as complex structure in 
momentum distributions. 

This work is supported by AFOSR (grant No. FA9550- 
07-1-0412) and FOCUS (NSF grant No. PHY-0114336). 
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